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Abstract 

We study the evolution of the ground state and the excitation spectrum of 
the two and three dimensional attractive (negative-?/) Hubbard model as the 
system evolves from a Cooper pair regime for U <C t, to a composite boson 
regime for U S> t. Our work is motivated by the observation that the high 
temperature superconductors, with their short coherence lengths and unusual 
normal state properties, may be in an intermediate coupling regime between 
these two limits. A mean field analysis of pairing, suitably generalized to 
account for a shift in the chemical potential, is known to be able to describe the 
ground state crossover as a function of U/t. We compute the collective mode 
spectrum using a generalized random phase approximation (RPA) analysis 
within the equations of motion formalism. We find a smooth evolution of the 
Anderson mode for weak coupling into the Bogoliubov sound mode for hard 
core bosons. We then include a long-range Coulomb interaction and show 
that it leads to a plasmon which again evolves smoothly from weak to strong 
coupling. 

PACS numbers: 74.20.Fg, 71.45. Gm 
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I. INTRODUCTION 



Since their discovery by Bednorz and MulleiS, six years ago, the cuprate superconductors 
have consistently shown appreciable deviations!!, both in their superconducting and normal 
state properties, from the conventional superconducting metals. One of the most remarkable 
characteristics of the high-Tc materials is their short coherence length which is a few times 
the lattice spacing so that the pairs of fermions are weakly overlapping in real space. This is 
in marked contrast with the BCS superconductors where the pair size is orders of magnitude 
larger then the lattice spacing, and the number of fermions within a Cooper pair is very 
large. 

The problem of the crossover from a BCS state with Cooper pairs to a condensate of 
composite bosons was first considered at T = by Leggetti and further analyzed at finite 
temperatures by Nozieres and Schmitt-Rink@. These ideas were taken up again in the context 
of the high T c superconductors where it was suggested! that these systems might properly 
be described as being in an intermediate regime between the BCS limit and the Bose limit. 
Very recently, the normal state in the intermediate coupling regime has been investigated!! 
by Monte Carlo simulations of the 2D attractive Hubbard model, with a view to studying 
deviations from Fermi liquid behavior. It was founJ that a degenerate Fermi system shows 
anomalous spin-correlations above T c , providing a natural qualitative explanation of the 
spin-gap behavior observedB in the Knight shift and NMR relaxation rate of several high-T c 
systems. 

In this paper we study the collective excitations in the attractive Hubbard model at 
T = as the ground state evolves from the BCS to the Bose regime. A study of the 
collective modes as a function of coupling is important for at least two reasons. The critical 
temperature in the BCS weak coupling regime, as well as the thermodynamic properties, 
are governed by the pair breaking excitations with an energy gap A which is exponentially 
small. However, in the strong coupling regime of composite bosons the energy gap becomes 
very large and pair breaking is no longer possible. T c is then controlled by the center of 
mass motion of the pairs, or in other words by the collective modes. In the intermediate 
regime, the thermodynamic properties will be controlled by a combination of pair breaking 
and collective modes. A study of the T = collective modes is then a first step towards 
understanding this complicated regime. 

A second reason for studying the collective modes is that, at least in the Bose regime, 
we know from Bogoliubov theory that their spectrum at long wavelengths depends in an 
essential way on the repulsive interaction between the bosons. While the crossover analyses, 
refered to above, were able to access the Bose regime while working with the constituent 
fermions, it is not apriori clear how the effective interactions between the pairs are accounted 
for in these treatments. 

We study the collective mode spectrum at T = using a generalized RPA formulation. 
We adapt the analysis of AndersonH and othersiil initially applied to weak coupling super- 
conductivity, and show that with some modifications it is capable of describing the evolution 
of the collective mode for all values of U / 1. While one might not expect the RPA to be valid 
in the strong coupling limit, we show that we recover the well known Bogoliubov resultlll 
for the sound velocity of a repulsive Bose gas in the U 3> t limit. Using the RPA as an 
interpolation scheme we find a smooth evolution of the collective mode spectrum in both 
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two and three dimensions. 

We extend the anaysis to charged superconductors in two and three dimensions. In 
weak coupling we find that the sound mode for neutral system is pushed up! to the plasma 
frequency in 3D; in 2D the plasmon has a y/q dispersion. We find that the plasmon evolves 
smoothly as a function of the attraction, and in the strong coupling, dense system limit we 
recover exactly the known plasma frequency for a dense charged Bose ga 

Some of the results of this paper have been presented without detailed derivation in a 
previous Rapid Communication^. Related results have been obtained by several authorfS0 
using different techniques. The detailed presentation given here is nevertheless of some in- 
terest, especially since we derive analytical results for the BCS and Bose limiting cases. 
Recently, the RPA has also been successfully applied to other crossover problems, for ex- 
ample, the evolution from itinerant to local moment antiferromagnetismlla, and excitonic 
collective modes in a Bose condensed electron-hole gas@. 



II. MEAN FIELD ANALYSIS 

Our starting point is the single band attractive Hubbard model on a <i-dimensional hyper- 
cubical lattice: we will focus on the 2D and 3D cases. The Hamiltonian 

H = -tJ2 4aCja n i\ n a - £ ( n *t + n a) (!) 

ijcr i i 

may be rewritten in momentum space as 

H = Y,( £ k - ^)4* c k* -JjY. c t+ g A'-Ql c k'l c ^ ( 2 ) 

k,a kk'q 

where Ek = — 2t X^=i cos(/c;a). The chemical potential /i will be adjusted to obtain the 
required band filling / = N/2M, where N is the average number of electrons and M the 
total number of sites. We will study the model at T = for arbitrary U/t > and 
< / < 1/2 where the ground state is expected to have superconducting off-diagonal long 
range order. We will not discuss the competition between superconductivity and CDW 
ordering^ at half-filling / = 1/2. 

We quickly review the mean field analysis, mainly to establish notation. The BCS 
"reduced" Hamiltonian H\> cs is that part of (|2|) which describes the interaction between 
pairs with zero center-of-mass momentum, ifbcs is thus obtained by retaining only the 
fcPhys.Rev.,ime _ pj ece Q f se cond term of (0). The BCS-Bogoliubov solution consists of 
determining the eigenoperators 7^ and 7fe CT of i^bcs- These operators define both the BCS 
ground state, via 7^ |$o) — 0, and the single-particle excitations 7^.|$o) with energy Ek- 

The first step is a linearization of the equations of motion for c^ a and c^ a with respect 
to |$o) where (c^ a Cka) and (cjkcj^i) are non- vanishing. We must retain the Hartree shift in 
the chemical potential 

fi = (i + fU (3) 

since it plays a crucial role for large U (see below). The second step is to diagonalize 
the linearized equations via the Bogoliubov transformation: = Uk-fw + v klhi anc ^ c -fci = 
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—Vklko + Uk^ki- Here u\ = l — v% — | (1 + £k/Ek), with £ fc = — /i, and the order parameter 
and the quasiparticle excitation energy are given by A = UY,k u kVk, and E k = (£| + A 2 ) 1 / 2 
respectively. 

Self-consistency is achieved by demanding that, for each value of the coupling U/t and 
filling /, A and the chemical potential \x satisfy the gap equation 

u v 1 m 

2M^[el + A 2 ]V2 ^ 
which is familiar from BCS theory, and the number equation 

/ = ^M?( 1_ [^ + A 2 ] 1/2 )' (5) 

The latter, which simply follows from 2 J2k v \ = is trivially solved in the BCS limit to 
obtain \i ~ ep, the noninteracting result. As the coupling grows, however, the occupation 
probabality in momentum space rife = 2v\ broadens significantly, corresponding to the for- 
mation of tightly bound pairs. As a result the chemical potential is strongly affectedi by 
the interactions (in addition to the trivial Hartree shift (^|). 

The gap and number equations may be solved analytically in the weak and strong cou- 
pling limitsi'S. (In the continuum limit in 2D an exact analytical solution is possible for 
arbitrary coupling si). For U/t <C 1 the chemical potential is at Ef, the pair size is much 
larger than the lattice spacing, as expected of a BCS ground state, and the energy gap 
has the usual essential singularity in U/t. In the opposite limit of U/t ^> 1, the pairs are 
on-site, the chemical potential is one-half the pair binding energy, and the ground state is a 
condensate of composite bosons. The expressions for jx and A in the strong coupling regime 
were initially given by Nozieres and Schmitt-Rinki; 

£ = ^(2/-l); A = UyffO. - f) (6) 

There is a smooth crossover between these rather different limits, as can be seen by a 
numerical solution of (|j) and ([5]). The solutions in 3D case are given in ref.Q and the 2D 
results are shown in Fig. la and b where we plot A and fi as a function of U for various 
fillings /. The gap to single-particle excitations is defined by 



-Egap = min^/ (e fc — /i) 2 + A 2 (7) 



where the minimum is to be found within the band, i.e., E& > —2dt. We thus find £ g „ = A 
provided £l lies within the band. However, once fl is below the bottom of the band™, we 

find E mv = ((-2dt + \fi\f + A 2 ) 17 '. 

To understand how a BCS-like analysis is able to describe a condensate of composite 
bosons for U 3> t, note that the (un-normalized) BCS ground state can be written as 

l^o) = (j2k V 9 fcCfc|C^ fc ^ |vac), with tpk = Uk/v k . Variationally, the mean field solution 
corresponds to an optimal choice of the internal pair wave-function (p, which in the strong 
coupling limit is an on-site singlet. The Pauli exclusion principle for the constituent fermions 
gives rise to a hard core repulsion for these composite bosons. As shown below, this will 
have an important effect on the collective mode spectrum. 
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III. GENERALIZED RPA 



At the mean field level we treated only pairs with zero center of mass momentum, by 
retaining only the ifbcs part of the full Hamiltonian H = H^, cs + H mt . In this section we 
shall treat within a generaized RPA the fluctuations introduced by H int which describes the 
interaction between the Bogoliubov quasiparticles. 

The RPA can be implemented in a variety of ways to study the collective excitations 
above the superconducting ground state. We use here the linearized equations of motion 
method which, in our view, is physically transparent and has a certain intuitive appeal. 
This method had been originally developed for the normal fermi liquid state by Bohm and 
Pinesil, and adapted to the superconducting state by Anderson!. 

To determine the collective mode spectrum we study the time evolution of density fluc- 
tuations Pfc CT (g) = cl +qa Cka- Its equation of motion is coupled to that of pairs with a finite 
center of mass momentum b k ^(q) = c k+q ^ct k ^ and b k ^(q) = c_ k ^ q ^c k ^, resulting from the 
particle-hole mixing due to the condensate. The first step in the generalized RPA is to 
linearize the equations of motion of these operators and the second to diagonalize them by 
finding the appropriate eigenoperators for the collective coordinates. 

Before proceeding with the calculation, it may be useful to write the terms appearing in 
the linearized equations of motion as diagrams. In Fig. 2 we show the first order processes 
involved in the full linearized equations of motion for c k+q ^ct k ^ (pair or p-p channel) and 
in Fig. 3 those for c k+ c ka (density fluctuation, or p-h channel). Note that no exchange 
scattering is allowed by the on-site (attractive) Hubbard interaction. In Fig. 2, (2a) is the 
usual particle-particle scattering vertex, and (2b) is just the Hartree self-energy term which 
renormalizes the chemical potential. The existence of the condensate (a non-zero expectation 
value of < b kcr (q = 0) >) mixes the p-p and p-h channels as shown in (2c) and (2d). 

We now plunge into the rather lengthy algebra of the equations of motion method. We 
use a derivation based on the very clear presentation of Bardasis and Schriefler0, adapting 
it to the lattice model, and retaining terms which allow us to work at arbitrary U. The weak 
coupling BCS analysis is considerably simplified by an approximate particle-hole symmetry 
arising due to the fact that only fermions in a thin shell which is symmetric about the Fermi 
surface are affected by the pairing. For arbitrary U one no longer has such a p-h symmetry 
which makes the calculation rather more complicated. 

We look at the time evolution of Jk+qalk,a, lk+q olkii anc ^ lk+ q ,i1k,o, (instead of working 
with bilinear products of c and c + ). Using the notation 

E P; q = E p+g + Ep (8) 

the Anderson- Rickayzen equations are given by the following commutation relations: 

k 

KP, O) Z! ?) (jfc+ g ,o7M - 7fc+ 9 ,i7fc,o 



2M 
U_ 
~2M 



■ n{p, q) n(k, q) (tm^otJi + 7fc+<?,i7fc,o) (9) 
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[H, 7 P + g ,i7 P ,o] = -^p, 9 7p+ 9 ,i7p,o + 2]T7 m (^> 9) Z m ( fc > ?) (7fc +g> o7ji + 7fc+ g ,i7fc,o) 



Z(P> 9) Z (7fe"+ 9j o7ji - 7fc+g,i7fc,o) 



2^ 



+ 2lf ^ ^ ^ ( 7 *+ff.° 7 M + 7fc +9.i7fc,o) (10) 
7fc +9 , CT 7fe, CT ] = (£fc+» - ^Jfc)7fc+ g , CT 7fc,<7- (11) 

We dropped all terms of the form 7 + 7 on the right hand side of the first two equations for 
reasons to be discussed shortly. The coherence factors are defined as 

l(k, q) = u k u k+q + v k v k+q 
m(k, q) = u k v k+q + v k u k+q 
n(k, q) = u k u k+q - v k v k+q 

p(k, q) = u k v k+q - v k u k+q (12) 

The next step is to diagonalize the above equations by finding eigenoperators which 
satisfy 

[H, fi + (q)} = Lj(q) fi + (q) (13) 

[H,n(q)] = -u>(q)ii(q). (14) 

The "renormalized" ground state |^o >, which differs from the mean field BCS ground state 
|<3?o > through the inclusion of zero-point collective excitations, is defined by 

li(q)\* o >=0. (15) 

In addition fi + (q) acting on > creates an excited state with an excitation energy u(q) > 
0. 

We see that ( [TT]) is already diagonal and the action of the eigenoperator on an initial 
state describes the scattering of a Bogoliubov quasiparticle already present in that state. 
There are no excitations present in the ground state, and lt+q alk,a\^ a >= 0. Since we will 
look at matrix elements of the equations of motion between > an d an excited state, 
this operator can be ignored in the subsequent discussion. As stated earlier, all the terms 



containing ^ k+ „ a ^k,a have already been dropped in @ and (|10|). 

The non-trivial eigenoperators n + (q) and fi(—q) must therefore be chosen as a linear 
combination of 7pt)_ 0j o7^i and 7 P +g,i7 P) o- Equivalently, we may write 



7fc+<2,o7jfc,i 



f(k,q)fM + (q)+f(k,qM-q) (16) 



7fc+g,l7fc,0 



g(k,q)n + (q) + g(k,q)fi(-q) 



(17) 



where the /'s and g's have to be determined. We substitute (|I6| ) and ([TT|) in the equations of 
motion @ and (|H]), and take matrix elements between the ground state > an d a state 
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\^(q) > with a single quantum of excitation with energy u(q). After some simple algebra 
we obtain 

M?) - E P ,q] fip, q) = m(p, q) Z(q) + - n(p, q) A(q) + ~ l(p, q) T(q) (18) 



H<?) + E P ,q] g(p, q) = ~ m(P, q) Z{q) - ^ n(p, q) A(q) + ^ Z(p, g) 

Note that / and g drop out using ( |l5l) , and ^(g), A(g) and are defined as 

U 



Z(q) 



A(q) 



2M 



lf( k ,q) + g(k,q)} 



u_ 

M 

u 



n ( k ^) lf( k ,q) + g(k,q)] 



r(?) = -^ [f(k,q)-g(k,q)] 



(19) 



(20) 
(21) 
(22) 



Using ( |l8l) and ([19]) we find 
1 



f + g 



uP{q) - El 



[2E Ptq m(p, q) Z(q) + E m n(p, q) A(q) + u{q) l(p, q) T(q)] (23) 



p.q 



f-g 



" 2 (q)-El q 



[2u{q) m(p, q) Z{q) + u(q) n(p, q) A{q) + E M l(p, q) T{q)} (24) 



Substituting these results in the expressions for Z(q), A(q) and T(q), we obtain the three 
coupled equations 



[i + u Ie,uAq)} Hq) + u L, n ,M) r(g) + 2Ui E>n , m {q) z{ q ) = o 

U I^ n j(q)A(q) + [1 + UI EM {q)\ T{q) + 2U I^ m {q) Z(q) = 
- I E ,m,n{q) A (<l) + \ Iu,,m,l(q) T(q) + [1 + U I E , m , m {<l)] %) = 



We have used (following ref.[0]) the notation 



1 ^ a(k,q)b(k,q)c(k,q) 
ia,b,c{q) - -Tj 2^ 



(25) 
(26) 

(27) 



(28) 



fe,q 



where a, 6, and c denote any one of the following quantities: the coherence factors defined 
12|) , the quasiparticle energy E^q defined in ®, or the excitation energy uj(q). 



m 



The collective mode spectrum u(q) is obtained by solving the secular equation 

[1 + U lE,n,n{q)] U L,n.l(q) 2U I E ,n,m(q) 

U I E , n , m (q) [1 + U I E ,u(q)} 2U I u j, m (q) 

(U/2) I E , m ,n(q) (U/2) L, m M I 1 + U I E ,mAl)} 



(29) 
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Before solving this equation analytically in several limiting cases, we can check a general 
feature of symmetry breaking. On general grounds we expect a Goldstone mode, i.e., a 
solution u(q) = corresponding to q = 0. (We shall show later that by including the 
long ranged Coulomb interaction, this mode will become massive via the Anderson- Higgs 
mechanism). By noticing from ( pp that as q — > and uj —>■ 0, I^.m^ and Iu,n,l vanish, while 
lE,n,n,lE,m,m, and lE,n,m remain finite, the secular equation reduces to 



l + UI Et u(0) = 1- 



U 1 







(30) 



This is identical to our gap equation 
of our results. 



), thus providing a nontrivial check of the consistency 



IV. WEAK COUPLING 



In weak coupling A is exponentially small and jl tends to the fermi energy e E . The 
integrals fl28|) are then peaked at £p. We use the particle-hole symmetry of the BCS limit 
to simplify the calculation. For q — 0, the products of coherence factors n(k, q) l(k, q) and 
n(k, q) m(k, q) are odd under change of sign of which leads to vanishing integrals for I Uin j 
and I E . n ,m- As a result we are left with the 2x2 determinant: 



l + UI EM {q) 
(U/2)I Utlim {q) 



2U I uJ t, m (q) 

1 + U I E ,m,m(q) 



o 



(31) 



The small q and to expansion of the various terms in ( |3TD is conveniently written in terms 
of four quantities: x = M" 1 ^^" 3 , V = M~* E k (V kt) 2 / E 3 k , w = M' 1 J2k^U/ E h and 
z = A 2 M _1 J2k(y 'kO 2 '/ ' and is given to leading order by 



l + UI EM {q) 



U 



d 



(3z + w — y)q 2 — xoj' 



(32) 



1 + U I E , m ,m{q) 



1--A 2 x 
2 



(33) 



u i u ,i, m (q) = (-^Ax)uj 

Substituing these results into (|3TD, we obtain the long wavelength dispersion relation 

'1 (1 - (U/2)A 2 x){3z + W-V) 1 1/2 



(34) 



OJ 



d 



X 



Q 



(35) 



where d is the dimensionality. 

The momentum sums are peaked about the Fermi level and may be estimated by integrals 
over a thin shell of thickness 2u c , such that A ^ u c ^ W = 2dt, centered around the Fermi 
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energy. We then find x = 2iV(0)/A 2 , y = 2iV v (0)/A 2 , w = -2a 2 N(0) log(2w c /A), and 
z = 4iV v (0)/3A 2 , where the density of states N(g) = (2vr)- d J d d k5(£ k - f) and 



N*(t) = / (05(^)^-0. (36) 

The only dependence on the (arbitrary) cutoff cj c is in u>, which, however, is negligible 
compared to the other terms in Q3"5| ) in the weak coupling limit where A is small. (In the 
continuum limit w is identically zero). 

The weak coupling, long wavelength dispersion is thus given by 

u{q) = [(v 2 F )/d] 1/2 (l-UN(0))V 2 q (37) 

where the mean squared Fermi velocity is defined by {vp) = N v (0)/N(0). 

The density, or filling, dependence of {vp) depends upon the band structure. For the 
tight binding model used here, we obtain 

{v 2 F ) _Jd d k{EUsm 2 {k l a))6{^l 



(2ta) 2 Jd<tk5(£ k ) 

and the dependence of {vp) on the filling / is plotted in Fig. 5. In general the /-dependence 
is smooth with the following exceptions. In the 3D case there is a sharp dip in the collective 
mode speed of sound at / = 0.125 due to a van Hove singularity. In the 2D case the speed 
of sound goes to zero as / — > 1/2 due to the nesting at half- filling. 

In the continuum limit, with a parabolic dispersion, J {vp) = pp/m and fl37|) reduces to 
Anderson's weak coupling resultli in 3D. The collective mode in weak coupling is essentially 
the same for the lattice and continuum models, as one might expect since the size of the 
bound pairs £o is much larger than the lattice spacing a. 



V. STRONG COUPLING 

In the strong coupling limit the fermions bind into on-site singlet pairs, and the low 
energy physics of the attractive Hubbard model can be obtained by a projection on to the 
non-singly occupied subspace. This exact mappingil results in a system of hard core bosons 
described by 

#bosc = £(^) (runj-b+bj). (39) 

The composite bosons move only via virtual ionization with an effective hopping amplitude 
% = —2t 2 /U. The hard core constraint is due to the Pauli principle for the constituent 
fermions, and in addition the bosons interact with a nearest neighbor repulsion V = 2t 2 /U. 
In the long wavelength limit, the collective excitation of a dilute {na^ 1) 3D Bose gas is 
the Bogoliubov sound moddll with dispersion u = (Ann-bas/m^) 1 ^ 2 q. The density of bosons 
is = f /a 3 and the effective mass is ~ — l/t^a 2 = U/2t 2 a 2 . In terms of the parameters 
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/, a, U and t which characterize the constituent fermions, the Bogoliubov dispersion is given 
by 



oo{q)^Jvf (U 2 a/U)q, 



(40) 



where we take the scattering length a s , describing the interactions between the bosons, to 
be of the order of the lattice spacing a. 

We now show that one obtains essentially the same result from strong coupling limit of 
the general expression (p9), which comes from an RPA analysis of the constituent fermions. 
Using the gap and the number equations, (f|) and we find A = U[f(l — f)] 1 ^ 2 [1 — da 2 } 
and fx = U(2f — 1) [1 + 2da 2 } /2, to leading order in a = 2t/U < 1. A similar expansion of 
the various quantities in (|29| ) (see appendix A for details of derivations), yields 



da' 



f- UI E , n , n = 4/(1 - /) [l + 4(1 - 6/ + Qf 2 )da' 
UI E , m , m = (2/ - l) 2 + 8/(1 - /) [l - 3(2/ - l) 2 
UI E ,n,m = (2/ - 1)[4/(1 - f)] 1 ' 2 [l - (5 - 6(2/ - 2) 2 )da< 

w(2/-i; 



U Iuj,n,l 



u 



We also find, independent of the filling, 

l + U I EU = a 2 q 2 a 2 /2 



l-2(2-3(2/-l) 2 )da 2 

-l + (l + 6(2/-l) 2 )rfa 2 



u 2 /U 2 



(41) 
(42) 
(43) 

(44) 
(45) 

(46) 



With the further simplification of the dilute limit / < 1, the various quantities become: 



1 + UI E , n , n = 4/[l + 4rfa 2 ], 1 + UIe, 



, m , m = 1 - lQfda 2 , UI E , n , m = -2fV 2 [l + da?], 
UI u , n ,i = -u[l + 2da 2 } /U, and U J w , m ,, = 2ui f^^ 2 [-1 + 7da 2 } /U. We then obtain the RPA 
result for the dilute (/ < 1), strong coupling (U/t ^> 1) limit 



(47) 



uj(q) = J4df(4t 2 a/U)q 



in (i-dimensions. This result is essentially the same as that expected for the dilute Bose gas 
(fiOf). The RPA gets both the square root density dependence and the inverse boson mass 
dependence of the speed of sound correctly. The only difference is an overall factor of the 
order of unity, which is not unexpected since we had simply used a s ~ a in (|40"D . Quite 
remarkably, starting with interacting fermions and using RPA we were able to reach the 
regime of hard core bosons in the strong coupling limit. 



VI. CROSSOVER 



Encouraged by the success of the RPA in strong coupling, we might expect it to be a 
reasonable interpolation scheme all the way from weak to strong coupling. For intermediate 
couplings we have numerically solved the RPA equation (^9) as a function of the coupling 
U/t, and for various fillings, / = 0.05,0.25,0.45. As input to these equations we have used 
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the numerical solutions for A and fx obtained from the mean field gap and number equations, 
(|) and @; (for the 2D case see Fig. 1). To evaluate the /c-sums in (|29"D we used slightly 
different procedures in 2D and 3D. In three dimensions we numerically evaluated the density 
of states N{£) and the function iV v (£) (see (|3~6|)) for the simple cubic nearest neighbor tight 
binding band structure. In 2D we used analytical expressions, in terms of complete elliptic 
integrals, for these density of states functions for the square lattice. (See appendix B for 
derivation of N v (£) in 2D). 

The collective mode velocity is plotted in Fig. 4a (for 2D) and Fig. 4b (for 3D) as a 
function of the coupling. These results supercede those given in Fig. 2 of our previous 
papeiH which contained a numerical error. The analytical results in the small U, obtained 
from (|37|) , are separated for clarity from the numerical curves in Fig. 4a and 4b. In strong 
coupling, the / = 0.05 result (|47|), is compared with the numerical result at the same 
filling. The / = 0.45 numerical result is compared in strong coupling with the analytical 
result for half filling / = 0.5 which is much easier to obtain. We see that at all fillings, 
the numerical results smoothly interpolate between the analytical weak and strong coupling 
results. The 1/U dependence of c is apparent in strong coupling, at all fillings. Note also 
the non-monotonic dependence on / in weak coupling as discussed above. 

VII. CHARGED SYSTEMS 

We now extend our analysis to account for the effects of the long ranged Coulomb in- 
teraction between fermions. It is well known that in the BCS limit the singular behavior 
of the Coulomb interaction in the long wavelength limit has a dramatic effect on the col- 
lective modes. Anderson first showed! that in three dimensions, within the weak coupling 
approximation the sound mode is pushed up to the plasma frequency which is well above 
the energy gap. This showed that the collective modes were in fact quite unimportant for 
charged superconductors, except to restore gauge invariance in BCS theory@. However as 
the coupling increases, we expect the energy gap to increase and the plasma frequency to 
decrease due to the increasing effective mass of the bound pairs in the strong coupling lattice 
model. In the limit of very large U the spectrum should tend to the plasma frequency of 
a charged bose gasB u = ^Ati n^e 2 /m-Q which goes to zero as U goes to infinity, since the 
mass of the boson m B ~ U/t 2 . 

The inclusion of the Coulomb interaction in the negative-?/ Hubbard model is necessarily 
somewhat ad-hoc. The Negative-f/ Hubbard Model is only an effective model where —U 
can be viewed as the screened interaction, and where the bare Coulomb interaction V c (q) is 
absent. In the RPA, screening is accomplished by attaching to each vertex a series of the 
polarization vertex shown in Fig. 3a. We see that the RPA equations of motion automatically 
screen the interaction vertex in the processes shown in Fig. 2c and 3a. Therefore we can 
use the bare interaction V c (q) in these two terms. However the interaction vertices in the 
remaining processes in Fig. 2 and 3 cannot be connected to the polarization vertex and thus 
are not screened within RPA. We must therefore use the screened interaction, namely — U, 
for these terms. This procedure has the effect of coupling the bare interaction V c to the 
particle-hole channel only, where V c is given by 
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Tr , N 2(d- l)vre 2 

KM = («) 

and a d is the unit cell volume in ci-dimensions. 

Incorporating these modifications in the equations of motion for and c^c, and trans- 
forming to the Bogoliubov representation, we find that the equations of motion (|29|) are 
modified in the following way 

[1 + UI E ^ n (q)} A(q) + UI u ^{q) T(q) + 2U I E , n , m (q) Z(q) = (49) 



U I E ,n, m (q)Hq) + [1 + UI E ,iM] r (?) + W I U)l>m (q) Z(q) = 



(50) 



- V c (q) I E ^ n {q)K{q) - V c {q) I u>m ,i(q) T(q) + [1 - 2V c {q) lE, m ,m(q)] Z(q) = (51) 

The reason why V c enters only in the last row of (p9| ), replacing the factors of U/2, can be 
traced to the collective coordinate Z of (|20|) which involves the coherence factor m = uv'+vu' 
and is thus related to the p-h channel. 

We will restrict ourselves here to the dense limit, i.e., close to half filling. The dilute 
regime is more complicated, since we know that, in the strong coupling limit, we would 
obtain a system of charged bosons, which at sufficiently low densities, would form a Wigner 
crystal. Thus one would expect a first order phase transition to a crystal with decreasing 
density, which lies beyond the scope of the present work. To simplify the algebra we will work 
at half-filling / = 0.5 assuming a state with off-diagonal long ranged order (and ignoring 
CDW ordering). Since \x = for / = 0.5 independent of U, there is an exact particle- hole 
symmetry for all U which greatly simplifies the calculation. Thus I Ujnt i = lE,n,m = an d we 
are left with the 2x2 determinant. 



1 + UIeM 2UI w , hm (q) 
-V c (q) Iu3,i,m(q) 1 - 2V c (q) I E ,m,m{q) 



The collective modes are obtained by solving the equation 



2V c (q) 



<-E,m,m 



u- 



T 2 

uj,Lm 







(52) 



1 + UI E ,IJL 

It is convenient to write I a ,b,c{q) = Iabc + ^ Ia,b,c{q), where J° fec = I a ,b,M = °)- Using 



(53) 



E,m,m 



1 v 2 A 2 

M V E k (u 2 - AEl) 



(54) 



o 1 V ^ A 

w ,i,m M ^ E k (u 2 - 4E 2 k ) 



(55) 



2 El, 



l E,l,l 



M^u 2 ~AEl 



(56) 
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and noticing that I®i m = (2A/Uu) (1 + U I^u), h is straightforward to show, using the 
gap equation (|), that 



U- 



(4° ^ 2 



,Lm) 



1 + UIlu 



0. 



(57) 



The dispersion relation then simplifies to 
1 = 2V c (q) 



4A 4A 2 



(58) 



We will study now the limiting case of weak and strong coupling in the 3D and 2D cases. 



A. 3D Weak coupling 

We expect here the collective mode to be pushed up to a plasma frequency much higher 
than the energy gap. The dominant contribution is given by 5 lE,m,m{q), which after some 
algebra can be written in the form 

SI E , m , m (q) - ^— J (59) 

Using the same tricks that were used in the weak coupling analysis of the neutral case 
(below equation (j35"D ), we see that the second term on the right hand side can be neglected 
for A — > and the first term gives the result 

5I E , m>m (q) = ^N v (0). (60) 



Substituting this into (|58"D we find 



UJ 



1 ~ ^e 2 (4)iV(0), (61) 



3 

which reduces to the plasma frequency u 2 = A7cne 2 /m for a parabolic dispersion. 



B. 3D Strong coupling 

In the strong coupling regime A 3> u, and the term (4A 2 /u 2 )8 Ie,i,i dominates eq. (|58"1) . 
After a small q expansion of 5 Ie,i,i the dispersion relation can be written as 

AV i V (W 
1 = k(9) ^m?-eT- (62) 

Using the explicit form of V c (q) we then obtain 
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Using the boson mass m B = U/2t 2 a 2 (as in the neutral case), charge ee = 2e, and density 
hb = l/2a 3 (since the constituent fermions are close to half filling), we can rewrite the above 
result as 



u 2 = 4 " " b6b (64) 

This is exactly the plasma frequency of a dense system of charged bosons0, thus providing 
another important check of the validity of RPA in the strong coupling regime at T = 0. 



C. Collective modes in 2D 

In two dimensions, V c (q) = 2ne 2 /a 2 q, and we expect the plasma frequency u p ~ ^Jq as 
q — > 0. To calculate the collective mode spectrum, one can then use a small u and small q 
expansion of (|53|). The calculation proceeds in a manner analogous to the neutral system 
calculation, and we find the dispersion relation 

-V c A 2 x(3z + w -y)q 2 = xu 2 (65) 



where we use the quantities x, y, z,w introduced below (pl|). 

In the BCS limit we use the expressions for x,y,w and z derived earlier for the weak 
coupling neutral case (see below |3~5|) ) and obtain the result 

uj 2 = 27re 2 (v 2 F )N(0)q (66) 

For comparison, it might be useful to look at the continuum results where N(0) = m/27r, 
(vp) = (pf/iti) 2 , and n = 2irp F so the (|66"D can be written as u = (27re 2 n/m) ^Jq. 



In the strong coupling Bose limit one can use either (p5|), or go back to ( pS[ ) with the 
additional simplification that uj <ti A. Both routes yield the same result 

U= { — ) q (67) 

This may be written in terms of the bose parameters m-Q = U/2t 2 a 2 , hb = l/2a 2 , and 

ee = 2e to obtain u = (2ne B n-B/mB) 1 ^ 2 y/q. 

We have numerically calculated the crossover in the 2D plasmon mode from the BCS 
regime to the Bose regime close to half-filling. The results are plotted in Fig. 6: we find a 
monotonically decreasing function of U which smoothly interpolates between the analytical 
results in weak (|66| ) and strong coupling flBTp. 
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VIII. CONCLUSION 



In conclusion we have studied the attractive Hubbard model in two and three dimensions 
which is the simplest lattice model that shows a crossover from weak coupling BCS theory 
to a strong coupling Bose condensed regime. Using a generalized RPA formulation for the 
collective excitations at T = we were able to analytically reproduce the known results in 
the BCS and the Bose limits, for both neutral and charged systems, and numerically show 
that there was a smooth evolution from one limit to the other. This suggests that the RPA 
provides a reasonable interpolation scheme in the intermediate regime. 

Several open questions remain, some of which are listed below. Generalization to layered 
superconductors should be of interest. Collective modes in layered systems in the weak 
coupling BCS regime have been studied recently^ and a formalism to study the crossover 
problem in layered superconductors has been developed by Cote and GrimnEJ. The general- 
ization to finite temperatures would also be of considerable interest. Recently Sa de Melo, 
Engelbrecht and one of the authors have used the functional integral method to address this 
problem in a continuum modelEI. The lattice model has some important differences from 
the continuum problem: in particular T c for the attractive Hubbard model is expected to 
grow in the BCS manner like exp(— t/U) for weak coupling but eventually drop likell t 2 /U 
in the Bose regime. It would be very interesting to derive a general expression for T c (U/t) 
and to determine its maximum value which necessarily occurs in the intermediate coupling 
regime. 
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I. APPENDIX A 



A. Part I 



In the first part of this appendix we shall calculate the quantities jl, A, and Eq = 
(jl 2 + A 2 ) 1 / 2 in strong coupling to order (t/U) 2 at arbitrary filling. We will often use the 
identities 



and 



k 



-L V el = 2dt 2 
M k 



(68) 



(69) 



We start by expanding the gap equation (|) to second order in Ek/E , which is at most 
of order t/U in the strong coupling limit: 



U v 1 

U 



2E n 



2E 2 ' 2E* 
I 3djl 2 w 2t^~ 



1 + d 3 



This can be solved for i?o to yield 



Now expanding the number equation 

2 ^ 1 + f-?t 



4 17' 



,2t, 



1 _ e\ - 2gfc/j 3el/i 2 



(70) 



(71) 



rearranging the equation and using the leading term of E expansion in the third term of 
the right hand side 



2/-1 



2/2 

77 



Again by inverting the equation, we find 

a = |(2/-i; 



1 - 2rf(|) 2 



l + 2rf(|) 2 



Substituing ( [73D in ( [70]) we find 



2t, 



l + rf(3(2/-l) 2 -l) (^ 



using A 2 = E 2 — ft 2 , and eqs.( [T3|, |T3J) we find 



A 



[4/(1 -/)] 1/2 (l-^(f) 2 ) 



(72) 



(73) 



(74) 



(75) 
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B. Part II 



Here we derive in detail the strong coupling expressions of the quantity 1 + U lE,n,n] the 
derivation of the other coefficients in the equations of motion follow is very similar. Since 
this quantity has a non- vanishing limit as q — > we can simply set q = at the outset, thus 
obtaining 



1 + U I E ,n,n = 1 



1 



u 



E 



2M k [(e k - jl) 2 + A 2 } 
U e\ + ft 2 - 2e k fi 



2M 



El 



3/2 

_ 3 e| - 2e fe /2 15gg./2 2 
2 £ 2 2E„ 4 



Using the property X^e*; = we find 

Up? 



1 + U I E ,n,; 



U 2 15/i 2 15/A 



•2 .\ //•.;; 2 .u/-;,, , 







2^ 4 



Using eqs( [73|, |74|) we get 

1 + 17 ii5,n,n = 1 - 



1 - 3 (3(2/ - I) 2 - l) d(|) a 



X 



(2/- l)^(L + 4d(|) 2 



4/(1 - /) 



2-15(2/ 
2t 



l + 4(l-6/ + 6/ 2 K-) 2 



Equation ( u&) is identical to eq.41 



15(2/-1) 4 U( 



2t, 



(76) 



II. APPENDIX B 



We will derive here the expression of the weighted density of states in 2-D defined by 

d 2 k 



NJe) 



(2vr) ; 



;(Ve k Y6(e-e k ) 



(77) 



This can be extracted from the advanced green's function!! using the relat 



ion 



N v {e) = -Im(Gt(e)) 

IT 



(78) 



G+(z) 



(2vr) 2 
(2ta) 2 
(2tt) 2 
{Waf 
(2vr) 2 



2 jv. fc ) 2 



d l k 



z- s k 



dk\ / dk 2 



sin 2 ki + sin 2 /c 2 

- 2£(cos k± + cos fe) 
cos 2 0i sin 2 02 

— W COS 01 COS 02 



(79) 
(80) 
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for a derivation of (j80|) see re E. We have now the product of two simple integrals which 
can be put in the form 



Gt(z) 



(Wa) 



>1 COS (pi 



W COS 0i 



(81) 



where W = it is half the bandwidth, a is the lattice spacing, and 

,2 JL 







sin 



2 



x — cos <p W cos (fti 

A contour integral evaluation of the integral above leads to 

-2tt 



x 



+ v 



X/ 



(82) 



(83) 



Substituing fl8"Bp into (FH]) we find 



—2; 



1-- / , 

71 JO 

1 - -E(W/z) 

71 



>yjl ~ (W/z) 2 COS 2 (j) 



(84) 



where E is the complete elliptic integral of the second kind. In the case where VF/kJ ^ 1 
(i.e e is taken within the band), one can use the analytic continuation of E given byc3 



E 



E(e) - iE(VT^) - e 2 K(e) + ie 2l K{s/l^) 



(85) 



where K is the complete elliptic integral of the first kind, and e is expressed in units of W. 
Taking the imaginary part Q85|) , substituing it into (|84"|) and using (|78|) we finally find 



NJe) = — 



7T Z 



E(VT zr £ f )-e ; 



K(VT 



(86) 



Equation fl8"B]) was used in the calculation of the energy gap A, the chemical potential /i, as 
well as the velocity of the collective mode c in 2D. For completeness we give also the density 
of states N(e) in 2D!. The figures of N(e) and N v (e) for both the 2D and 3D cases are 
shown in figure 7. 



N(e) 



(87) 
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FIGURES 



FIG. 1. (a) The order parameter A and b) the chemical potential fi, as functions of the 
coupling U for 2D. All energies are plotted in units of the half-bandwidth W = 4t.The dashed, 
solid and dotted lines correspond to filling factors / = 0.45, 0.25, and 0.05 respectively. 

FIG. 2. Diagrammatic representation of various terms arising to first order in the interaction 
in the pair, or particle-particle channel. The dashed line represents the interaction U, and the 
black circle represents the condensate. 

FIG. 3. Diagrammatic representation of various terms arising to first order in the interaction 
in the density fluctuation, or particle-hole channel. The dashed line represents the interaction U, 
and the black circle represents the condensate. 

FIG. 4. Collective Mode velocity as a function of coupling U for various fillings in (a) two 
and (b) three dimensions. W = 2dt is half the bandwidth, d the dimensionality, and a the lattice 
spacing. The dashed, solid and dotted lines correspond to / = 0.45, 0.25, and 0.05 respectively. 
The results of the numerical evaluation are separated by breaks in the curve from the weak and 
strong coupling results obtained anaytically. 

FIG. 5. Fermi velocity in 2D and 3D as a function of filling in the weak coupling limit. The 
sharp feature near / = 0.21 for the 3D case is a result of a van Hove singularity. 

FIG. 6. Collective Mode velocity as a function of coupling U for the 2D charged system in the 
dense limit (/ ~ 0.5). The result of the numerical evaluation is separated by breaks in the curve 
from the weak and strong coupling results obtained anaytically." 

FIG. 7. The density-of-states functions N(e) and N v (e) for (a) 2D and (b) 3D tight binding 
model. 
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